# Analysis file for:
# Carnegie, Clark, and Zucker, "Global Governance under Populism"
# World Politics

# Software: R 4.3.2 / RStudio 2024.04.1+748 / macOS Sonoma 14.5
# Hardware: 2019 MacBook Pro, 2.4 GHz 8-Core Intel Core i9, 64 GB

# Loading packages
pacman::p_load(tidyverse,
               broom,
               data.table,
               interflex,
               janitor,
               lfe,
               readxl)

# ------------------ #
# -> DATA LOADING ####
# ------------------ #

WDI_MISS <- read_csv("WDI_MISS.csv")
GHG_GAP1 <- read_csv("GHG_GAP_annex1.csv")
WDI_BACKFILL <- read_csv("wdi_stack_valid_agg_by_var_year.csv")
WDI_BACKFILL_BY_POPULISM <- read_csv("backfilling_by_populism.csv")
WDI_BACKFILL_BY_POPULISM_SCI <- read_csv("backfilling_by_populism_sci.csv")

# --------------- #
# -> MAIN TEXT ####
# --------------- #

# Table 1(1): Scientific Missingness ~ Populism
WDI_MISS %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l1 | country_code + year | 0 | country_code,
       data = .) -> reg_sci_all00
tidy(reg_sci_all00)

# Table 1(2): Scientific Missingness ~ Populism
WDI_MISS %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l1 + polity2 | country_code + year | 0 | country_code,
       data = .) -> reg_sci_all01
tidy(reg_sci_all01)

# Table 1(3): Scientific Missingness ~ Populism
WDI_MISS %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l1  +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_sci_all02
tidy(reg_sci_all02)

# Table 2(1): Scientific Missingness (Raw vs. Estimated/Imputed) ~ Populism
WDI_MISS %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l1 | country_code + year | 0 | country_code,
       data = .) -> reg_sci_raw00
tidy(reg_sci_raw00)

# Table 2(2): Scientific Missingness (Raw vs. Estimated/Imputed) ~ Populism
WDI_MISS %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l1 + polity2 | country_code + year | 0 | country_code,
       data = .) -> reg_sci_raw01
tidy(reg_sci_raw01)

# Table 2(3): Scientific Missingness (Raw vs. Estimated/Imputed) ~ Populism
WDI_MISS %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l1  +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_sci_raw02
tidy(reg_sci_raw02)

# Table 2(4): Scientific Missingness (Raw vs. Estimated/Imputed) ~ Populism
WDI_MISS %>%
  felm(SCIENCE_NA_IMP_sd ~ fpopulism_aer_l1  +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_sci_imp02
tidy(reg_sci_imp02)

# Table 3(1): Emissions Gap ~ Populism
GHG_GAP1 %>%
  felm(diff ~ fpopulism_aer_l1 | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_annex0
tidy(reg_ghg_annex0)

# Table 3(2): Emissions Gap ~ Populism
GHG_GAP1 %>%
  felm(diff ~ fpopulism_aer_l1 + polity2 | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_annex1
tidy(reg_ghg_annex1)

# Table 3(3): Emissions Gap ~ Populism
GHG_GAP1 %>%
  filter(execrlc != -999) %>% #removing NAs
  felm(diff ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc == 1) +
         log1p(gdp_per_capita) +
         IMF +
         ff.energy.consumption +
         ag.va | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_annex2
tidy(reg_ghg_annex2)

# Figure 1: GHG Reporting Gap ~ WDI Missingness

## Left-hand plot
GHG_GAP1 %>%
  mutate(ghg_gap_sd = scale(abs(diff)),
         wdi_miss_sd = scale(SCIENCE_NA)) %>%
  ggplot(aes(x = wdi_miss_sd, y = ghg_gap_sd)) +
  geom_point(aes(color = as.character(fpopulism_aer_l1),
                 alpha = as.character(fpopulism_aer_l1))) +
  scale_color_discrete("Populist (t-1)") +
  scale_alpha_discrete("Populist (t-1)") +
  xlim(-2, 1.2) +
  theme_classic() +
  xlab("All WDI Missingness (Standardized)") +
  ylab("GHG Reporting Gap (Standardized)") +
  theme(axis.title.x = element_text(size = 10),
        axis.title.y = element_text(size = 10))

## Right-hand plot
GHG_GAP1 %>%
  mutate(ghg_gap_sd = scale(abs(diff)),
         wdi_miss_sd = scale(SCIENCE_NA_RAW)) %>%
  ggplot(aes(x = wdi_miss_sd, y = ghg_gap_sd)) +
  geom_point(aes(color = as.character(fpopulism_aer_l1),
                 alpha = as.character(fpopulism_aer_l1))) +
  scale_color_discrete("Populist (t-1)") +
  scale_alpha_discrete("Populist (t-1)") +
  xlim(-2, 1.2) +
  theme_classic() +
  xlab("Raw WDI Missingness (Standardized)") +
  ylab("GHG Reporting Gap (Standardized)") +
  theme(axis.title.x = element_text(size = 10),
        axis.title.y = element_text(size = 10))

# Figure 2: Missingness Rate ~ Years Since Variable Year

WDI_BACKFILL %>%
  filter(year >= 2004) %>%
  group_by(years_till_version) %>%
  summarize(na_rate = mean(na_rate, na.rm = TRUE)) %>%
  ggplot(aes(x = years_till_version, y = na_rate)) +
  geom_line() +
  theme_classic() +
  ylim(.2, 1) +
  xlab("Years Since Variable Year") +
  ylab("Missingness Rate (Mean)")

# Figure 3: Missingness Rate ~ Years Since Variable Year (By Populist / Non-Populist Government)

## Left-hand plot
WDI_BACKFILL_BY_POPULISM %>%
  ggplot(aes(x = years_since, y = na_rate, group = fpopulism_aer_l1, lty = as.character(fpopulism_aer_l1))) +
  geom_line() +
  scale_linetype_discrete("Populist Government in Variable Year") +
  theme_classic() +
  ggtitle("All Variables") +
  ylab("Missingness Rate") +
  xlab("Years Since Variable Year") +
  ylim(.2, 1) +
  theme(legend.position = "bottom")

## Right-hand plot
WDI_BACKFILL_BY_POPULISM_SCI %>%
  filter(sci == TRUE) %>%
  ggplot(aes(x = years_since, y = na_rate, group = fpopulism_aer_l1, lty = as.character(fpopulism_aer_l1))) +
  geom_line() +
  scale_linetype_discrete("Populist Government in Variable Year") +
  theme_classic() +
  ylab("Missingness Rate") +
  xlab("Years Since Variable Year") +
  ylim(.2, 1) +
  theme(legend.position = "bottom") +
  ggtitle("Scientific Variables")

# ---------------- #
# -> APPENDICES ####
# ---------------- #

# Table D7: Weighted Missingness

WDI_MISS %>%
  felm(SCIENCE_NA_WEIGHT_sd ~ fpopulism_aer_l1  +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_weight1
tidy(reg_weight1)

WDI_MISS %>%
  felm(SCIENCE_NA_RAW_WEIGHT_sd ~ fpopulism_aer_l1  +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_weight2
tidy(reg_weight2)

WDI_MISS %>%
  felm(SCIENCE_NA_IMP_WEIGHT_sd ~ fpopulism_aer_l1  +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_weight3
tidy(reg_weight3)

# Table E8: Two Years Pre/Post Populist Entry

wdi_beforeAfter <- WDI_MISS %>%
  group_by(country_code) %>%
  mutate(populism_tenure = replace_na(populism_tenure, 0),
         post_pop = as.integer(dplyr::lag(populism_tenure, 2) > 0 & dplyr::lag(populism_tenure, 1) == 0 & populism_tenure == 0),
         pre_pop = as.integer(dplyr::lead(populism_tenure, 2) == 1 & populism_tenure == 0), # 2016
         first_pop = as.integer(populism_tenure == 2),
         last_pop = as.integer(dplyr::lead(populism_tenure, 1) > 0 & dplyr::lead(populism_tenure, 2) == 0)) %>%
  ungroup()

wdi_beforeAfter %>%
  filter(pre_pop == 1 | first_pop == 1) %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l2 | country_code | 0 | country_code,
       data = .) -> prepost1
tidy(prepost1)

wdi_beforeAfter %>%
  filter(pre_pop == 1 | first_pop == 1) %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l2 | country_code | 0 | country_code,
       data = .) -> prepost2
tidy(prepost2)

wdi_beforeAfter %>%
  filter(post_pop == 1 | last_pop == 1) %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l2 | country_code | 0 | country_code,
       data = .) -> prepost3
tidy(prepost3)

wdi_beforeAfter %>%
  filter(post_pop == 1 | last_pop == 1) %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l2 | country_code | 0 | country_code,
       data = .) -> prepost4
tidy(prepost4)

# Table E9: Three-Year Lag

WDI_MISS %>%
  filter(fpopulism_aer_l1 == 0) %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l3 +
         polity2_l3 +
         rightwing_l3 +
         gdp_per_capita_ln_l3 +
         IMF_l3 | country_code + year | 0 | country_code,
       data = .) -> reg_l3_1
tidy(reg_l3_1)

WDI_MISS %>%
  filter(fpopulism_aer_l1 == 0) %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l3 +
         polity2_l3 +
         rightwing_l3 +
         gdp_per_capita_ln_l3 +
         IMF_l3 | country_code + year | 0 | country_code,
       data = .) -> reg_l3_2
tidy(reg_l3_2)

WDI_MISS %>%
  filter(fpopulism_aer_l1 == 0) %>%
  felm(SCIENCE_NA_IMP_sd ~ fpopulism_aer_l3 +
         polity2_l3 +
         rightwing_l3 +
         gdp_per_capita_ln_l3 +
         IMF_l3 | country_code + year | 0 | country_code,
       data = .) -> reg_l3_3
tidy(reg_l3_3)

# Table F10: Removing Outliers (WDI Missingness)

WDI_MISS %>%
  filter(SCIENCE_NA_sd >= -2 &
           SCIENCE_NA_sd <= 2) %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         gdp_per_capita_ln +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_outlier1
tidy(reg_outlier1)

WDI_MISS %>%
  filter(SCIENCE_NA_RAW_sd >= -2 &
           SCIENCE_NA_RAW_sd <= 2) %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         gdp_per_capita_ln +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_outlier2
tidy(reg_outlier2)

WDI_MISS %>%
  filter(SCIENCE_NA_IMP_sd >= -2 &
           SCIENCE_NA_IMP_sd <= 2) %>%
  felm(SCIENCE_NA_IMP_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         gdp_per_capita_ln +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_outlier3
tidy(reg_outlier3)

# Table F11: Removing Outliers (Emissions Gap)

.ghg_quality_ub <- mean(GHG_GAP1$diff, na.rm = TRUE) + 2*sd(GHG_GAP1$diff, na.rm = TRUE) # +2SD bound
.ghg_quality_lb <- mean(GHG_GAP1$diff, na.rm = TRUE) - 2*sd(GHG_GAP1$diff, na.rm = TRUE) # -2SD bound

GHG_GAP1 %>%
  filter(execrlc != -999) %>% # removing NAs
  filter(diff <= .ghg_quality_ub & diff >= .ghg_quality_lb) %>%
  felm(diff ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc == 1) +
         log1p(gdp_per_capita) +
         IMF +
         ff.energy.consumption +
         ag.va | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_outlier
tidy(reg_ghg_outlier)

# Table G12: Non-Scientific Missingness

WDI_MISS %>%
  felm(NONSCIENCE_NA_sd ~ fpopulism_aer_l1 | country_code + year | 0 | country_code,
       data = .) -> reg_nonsci1
tidy(reg_nonsci1)

WDI_MISS %>%
  felm(NONSCIENCE_NA_sd ~ fpopulism_aer_l1 + polity2 | country_code + year | 0 | country_code,
       data = .) -> reg_nonsci2
tidy(reg_nonsci2)

WDI_MISS %>%
  felm(NONSCIENCE_NA_sd ~ fpopulism_aer_l1  +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_nonsci3
tidy(reg_nonsci3)

# Table H13: Dropping U.S. (WDI Missingness)

WDI_MISS %>%
  filter(country_code != "USA") %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_dropusa1
tidy(reg_dropusa1)

WDI_MISS %>%
  filter(country_code != "USA") %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_dropusa2
tidy(reg_dropusa2)

WDI_MISS %>%
  filter(country_code != "USA") %>%
  felm(SCIENCE_NA_IMP_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_dropusa3
tidy(reg_dropusa3)

# Table H14: Dropping U.S. (Emissions Gap)

GHG_GAP1 %>%
  filter(execrlc != -999) %>% # removing NAs
  filter(country_code != "USA") %>%
  felm(diff ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc == 1) +
         log1p(gdp_per_capita) +
         IMF +
         ff.energy.consumption +
         ag.va | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_dropusa
tidy(reg_ghg_dropusa)

# Table I15: Dropping Post-2015 (WDI Missingness)

WDI_MISS %>%
  filter(year <= 2015) %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_drop2015post1
tidy(reg_drop2015post1)

WDI_MISS %>%
  filter(year <= 2015) %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_drop2015post2
tidy(reg_drop2015post2)

WDI_MISS %>%
  filter(year <= 2015) %>%
  felm(SCIENCE_NA_IMP_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_drop2015post3
tidy(reg_drop2015post3)

# Table I16: Dropping Post-2015 (Emissions Gap)

GHG_GAP1 %>%
  filter(execrlc != -999) %>% # removing NAs
  filter(year <= 2015) %>%
  felm(diff ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc == 1) +
         log1p(gdp_per_capita) +
         IMF +
         ff.energy.consumption +
         ag.va | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_drop2015post
tidy(reg_ghg_drop2015post)

# Table J17: Blair Institute Measure (WDI Missingness)

WDI_MISS %>%
  felm(SCIENCE_NA_sd ~ populism  +
         polity2  | country_code + year | 0 | country_code,
       data = .) -> reg_blair1
tidy(reg_blair1)

WDI_MISS %>%
  felm(SCIENCE_NA_sd ~ populism  +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_blair2
tidy(reg_blair2)

# Table J18: Blair Institute Measure (Emissions Gap)

GHG_GAP1 %>%
  filter(execrlc != -999) %>% # removing NAs
  felm(diff ~ populism +
         polity2 +
         I(execrlc == 1) +
         log1p(gdp_per_capita) +
         IMF +
         ff.energy.consumption +
         ag.va | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_blair
tidy(reg_ghg_blair)

# Table K19: Additional Covariates (WDI Missingness)

WDI_MISS %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         fossil_fuels +
         agriculture_va +
         net_oda +
         gdp_growth_annual_pct + 
         unemployment_ilo_estimate +
         crisis | country_code + year | 0 | country_code,
       data = .) -> reg_additional_covariate1
tidy(reg_additional_covariate1)

WDI_MISS %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         fossil_fuels +
         agriculture_va +
         net_oda +
         gdp_growth_annual_pct + 
         unemployment_ilo_estimate +
         crisis | country_code + year | 0 | country_code,
       data = .) -> reg_additional_covariate2
tidy(reg_additional_covariate2)

WDI_MISS %>%
  felm(SCIENCE_NA_IMP_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         fossil_fuels +
         agriculture_va +
         net_oda +
         gdp_growth_annual_pct + 
         unemployment_ilo_estimate +
         crisis | country_code + year | 0 | country_code,
       data = .) -> reg_additional_covariate3
tidy(reg_additional_covariate3)

# Table K20: Length of Populist Spell

WDI_MISS %>%
  mutate(populism_tenure = replace_na(populism_tenure, 0)) %>% # coding non-populist country-years as 0-year spell lengths
  felm(SCIENCE_NA_sd ~ populism_tenure +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         fossil_fuels +
         agriculture_va +
         net_oda +
         gdp_growth_annual_pct + 
         unemployment_ilo_estimate +
         crisis | country_code + year | 0 | country_code,
       data = .) -> reg_spell1
tidy(reg_spell1)

WDI_MISS %>%
  mutate(populism_tenure = replace_na(populism_tenure, 0)) %>% # coding non-populist country-years as 0-year spell lengths
  felm(SCIENCE_NA_RAW_sd ~ populism_tenure +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         fossil_fuels +
         agriculture_va +
         net_oda +
         gdp_growth_annual_pct + 
         unemployment_ilo_estimate +
         crisis | country_code + year | 0 | country_code,
       data = .) -> reg_spell2
tidy(reg_spell2)

WDI_MISS %>%
  mutate(populism_tenure = replace_na(populism_tenure, 0)) %>% # coding non-populist country-years as 0-year spell lengths
  felm(SCIENCE_NA_IMP_sd ~ populism_tenure +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         fossil_fuels +
         agriculture_va +
         net_oda +
         gdp_growth_annual_pct + 
         unemployment_ilo_estimate +
         crisis | country_code + year | 0 | country_code,
       data = .) -> reg_spell3
tidy(reg_spell3)

# Table L21: Nationalism Measure

WDI_MISS %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         natcon | country_code + year | 0 | country_code,
       data = .) -> reg_nationalism1
tidy(reg_nationalism1)

WDI_MISS %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         natcon | country_code + year | 0 | country_code,
       data = .) -> reg_nationalism2
tidy(reg_nationalism2)

WDI_MISS %>%
  felm(SCIENCE_NA_IMP_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         natcon | country_code + year | 0 | country_code,
       data = .) -> reg_nationalism3
tidy(reg_nationalism3)

# Table M22: World Bank Conditionalities

WDI_MISS %>%
  felm(SCIENCE_NA_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         conditions | country_code + year | 0 | country_code,
       data = .) -> reg_condition1
tidy(reg_condition1)

WDI_MISS %>%
  felm(SCIENCE_NA_RAW_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         conditions | country_code + year | 0 | country_code,
       data = .) -> reg_condition2
tidy(reg_condition2)

WDI_MISS %>%
  felm(SCIENCE_NA_IMP_sd ~ fpopulism_aer_l1 +
         polity2 +
         I(execrlc2 == 1) +
         log1p(gdp_per_capita) +
         IMF +
         conditions | country_code + year | 0 | country_code,
       data = .) -> reg_condition3
tidy(reg_condition3)

# Figure N5: Populism x Polity Interaction Plot

interflex(data = as.data.table(WDI_MISS),
          Y = "SCIENCE_NA_sd",
          X = "polity2",
          D = "fpopulism_aer_l1",
          Z = c("gdp_per_capita_ln", "IMF", "rightwing"),
          FE = c("country_code", "year"),
          cl = "country_code",
          estimator = "binning",
          full.moderate = F,
          nbins = 3,
          na.rm = TRUE,
          Ylabel = "Science WDI Missingness",
          Dlabel = "Populism",
          Xlabel = "Polity2",
          theme.bw = TRUE) -> int01

int01$figure

# Table N23: Populism x Polity2 Interaction Table

int01$est.bin

# Figure N6: Populism x GDP Per Capita Interaction Plot

interflex(data = as.data.table(WDI_MISS),
          Y = "SCIENCE_NA_sd",
          X = "gdp_per_capita_ln",
          D = "fpopulism_aer_l1",
          Z = c("polity2", "IMF", "rightwing"),
          FE = c("country_code", "year"),
          cl = "country_code",
          estimator = "binning",
          full.moderate = F,
          nbins = 2,
          na.rm = TRUE,
          Ylabel = "Science WDI Missingness",
          Dlabel = "Populism",
          Xlabel = "GDP per capita (ln)",
          theme.bw = TRUE) -> int02

int02$figure

# Table N24: Populism x GDP Per Capita Interaction Table

int02$est.bin

# Table O25: Emissions Gap ~ WDI Missingness

GHG_GAP1 %>%
  mutate(ghg_gap_sd = scale(abs(diff)),
         wdi_miss_sd = scale(SCIENCE_NA)) %>%
  felm(ghg_gap_sd ~ wdi_miss_sd*fpopulism_aer_l1 | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_wdi1
tidy(reg_ghg_wdi1)

GHG_GAP1 %>%
  mutate(ghg_gap_sd = scale(abs(diff)),
         wdi_miss_sd = scale(SCIENCE_NA)) %>%
  filter(execrlc != -999) %>% # removing NAs
  felm(ghg_gap_sd ~ wdi_miss_sd*fpopulism_aer_l1 + 
         polity2 +
         I(execrlc == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_wdi2
tidy(reg_ghg_wdi2)

GHG_GAP1 %>%
  mutate(ghg_gap_sd = scale(abs(diff)),
         wdi_miss_sd = scale(SCIENCE_NA_RAW)) %>%
  felm(ghg_gap_sd ~ wdi_miss_sd*fpopulism_aer_l1 | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_wdi3
tidy(reg_ghg_wdi3)

GHG_GAP1 %>%
  mutate(ghg_gap_sd = scale(abs(diff)),
         wdi_miss_sd = scale(SCIENCE_NA_RAW)) %>%
  filter(execrlc != -999) %>% # removing NAs
  felm(ghg_gap_sd ~ wdi_miss_sd*fpopulism_aer_l1 + 
         polity2 +
         I(execrlc == 1) +
         log1p(gdp_per_capita) +
         IMF | country_code + year | 0 | country_code,
       data = .) -> reg_ghg_wdi4
tidy(reg_ghg_wdi4)

